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ABSTRACT 


Motion control of underway replenishment operations is 
achieved through the use of sliding mode control with a Linear 
Quadratic Gaussian compensator design. External disturbances 
include first-order wave force and moment as well as slowly 
varying interaction forces and moments between the two ships. 
Feedback control is used to provide adequate stability of 
motions while feedforward control with disturbance estimation 
and compensation achieves the desired steady state accuracy. 
The results demonstrate that satisfactory path keeping during 
operations can be maintained for various ship proximity 


distances and environmental conditions. 
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I. INTRODUCTION 


A. AIM OF THIS STUDY 

Accurate course keeping is one of the most important 
problems for ship maneuvering and especially in an underway 
replenishment in an cpen sea. There is an appare it danger in 
such an operation due to the complex effects of the seaway, 


and the conventionally slow control loop indicated in Figure 


1.1. 


Rudder angle Indication 
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Figure 1.1 Close Loop Control System for Ship Maneuvering 


In particular, during underway replenishment the ship is 
subject to a variety of forces and environmental disturbances. 


The most important of these are the interaction forces and 








moments that are generated as a result of the two ships in 


proximity with one another [Ref.7] [Ref.8]. These forces depend 
on the relative positions of the two ships and both their 
Magnitude and sign is not easily predictable by analytical 
techniques In a random seaway there exist also first-order 
irregular wave forces and second-order wave drift forces. The 
latter appear in the form of long term s’owly varying 
excitations and do not have a significant effect on lateral 
separation control[Re:.8]. In this thesis we focus our 
attention on the effects of the interaction forces and first- 
order wave forces on the system response. We employ a sliding 
mode feedback control law due to its robustness 
characteristics with regard to unmodeled dynamics and 
disturbances. The sliding surface is based on the Linear 
Quadratic Regulator. Partial staze measurement is assumed and 
a Kalman filter is designed to provide an accurate estimate of 
the unmeasurable states and to minimize the effects of 
measurement noises. Disturbance estimation and compensation is 
used in order to ensure the desired steady state accuracy. A 
schematic block diagram of the control design is shown in 


Figure 1.2. 
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fagueis: 2.2 Automatic Contro) of Ship Replenishment at Sea 

















B. THESIS OUTLINE 
In Chapter II, the LOR method of sliding mode controller 
design is developed and simulated at calm sea. Chapter III is 


conserned with the Kalman filter design. Simulations of the 


UNREP with open sea effects are presented in Chapter III. 








II. MATHEMATICAL MODEL IN THE HORIZONTAL PLANE 


A. NONLINEAR STEERING EQUATIONS IN THE OPEN SEA 

The UNREP operation at sea is schematically depicted in 
Figure 2.1 where Y and N are the interaction forces and 
moments. The leading ship keeps its course unchanged and the 
tracking ship keeps the separation distance constant, normally 
less then 100 feet (depending on the type of ship). 

Figure 2.6 and Figure 2.7 show the interaction forces and 
moments on the tracking ship due to the interaction effect of 
the leading ship. The interaction force Y(A,B) acts through 
the ship’s center of gravity and the moment N(A,B) about the 
center of gravity in the horizontal plane. The two curves used 
in this study are for a typical Mariner ship at a nominal 
speed of 15 knots and for different relative positions. In 
applying these curves to the leading ship, the interaction 
force and moment need to be changed to -Y(-A,B) and -N(-A,B), 
respectively. For each figure, the curves for B=50 ft and 
B=100 ft were determined from experimental model testing 
results[Ref.6] [Ref.8], and the curves in between were 
determined by interpolation. The data were also extrapolated 
up to B=150 ft. A and B are continually read by the simulation 


program as input data to a two dimensional interpolation 











routine which calculates the base line interaction forces and 
moments acting on the ship. 

The coordinate system used in this study follows the 
standard system from Principles of Naval Architecture [Ref.11]. 
In practice, the leading ship keeps constant heading, while 
the tracking ship maintains a constant lateral separation 
distance. Different configurations and the directions of 
variations of the interaction force and moment as the tracking 
ship moves from -524 ft to 524 ft are shown in Figure 2.2. 

The primary consern for UNREP simulation is related to the 
Space coordinates (longitudinal and lateral separation 
distance, yaw angle, and yaw rate). So the equations of motion 
are first expressed in ship velocity coordinates (u, v) and 
then transformed to space coordinates(x, y) (Figure 2.3). The 


transformation equations are 


xX=ucosy-vsiny (2.1) 


y=usinyp+vcosy (2.2) 


The ship’s coordinates are assumed to be at the ship's 


center of gravity (x,=0, yg=0 ). 









SHIP 


| TRACKING 


LEADING SHIP 


.1 Relative Fosition of 
Passes the Leading Sh 


Pp 
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* Note : 
A: longitudinal separation distance. 
BE: lateral separation distance. 
Arrows indicate positive directions according to 


established sion convention [Ret.i2). 


~~) 





Figure 2.2 Relative Positions of The Two Ships for A=-524 
ft, 0, and +524 ft and Corresponding Interaction 
Forces and Moments. 


* Note: LS: Leading Ship. 


Tes 2 2aCkiag Shea. 





HiVd SdiHs 
SNOINVINVISNI 
OL INJONVI 
YOLIIA ALIDONIA ete 





‘JIVNIGUOO) 

-NOULISOd 

WNIQNLIONOI 
N 


vl 





ee 








7 JIVNIQUOO) 
NOILISOd 
ISHIASNVS! 






: Hivd \ 
S.dlHS \ 
\ 














o SIXV JDVdS 0 









Sr WIL lV diHS 4O 
ALIAVSS 4O Y31N3) JO NOILUSOYU 








B. MANEUVERING MATHEMATICAL MODEL 

The nonlinear maneuvering equations of motion and the 
hydrodynamic interaction forces and moments that act on the 
ship’s hull are modeled for 15 knots nondimensional forward 
speed. In an UNREP control simulation, the ability to manually 
or automatically control the rudder and propeller speed must 
be present. Only rudder input is considered here, and the ship 
forward speed is assumed constant by the propulsion control. 
The basic mariner study ship characteristics are presented at 


Table 2.1. 


Table 2.1 Characteristics of Mariner-type Study Ship 











Length between perpendicular 

Beam 

Draft 
(15120 t) 

[aie eee! 


Block coefficient (C,) 


* The ship’s coordinate system is assumed to be at the 
ship's center of gravity(xg, yg= 0). 
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The goal of this study is to present an alternative way of 


steering control based on sliding modes. For this reason the 
Simulation does not concentrate on the derivation of ship 
dynamic and hydrodynamic effects, but on the results of 
control and observer feedback. Only the tracking ship is the 
major consideration, and we assume a reliable lateral distance 
measurement. The ship hydrodynamic maneuvering model for the 
Study ship consists of the nonlinear equations of motion in 
the horizontal plane (sway and yaw). The hydrodynamic 
coefficients are presented in Table 2.3, and the 


nondimensional equations of motion are: 


SWAY EQUATION: 


(m-Y,) v+ (mx,-Y,) L=F,(u,v,r,6) 2 3) 


YAW EQUATION: 


(mx,-N,) +(I,-N,) f=F,(u,v,r,6) (2.4) 


where 


Bae, 2B) SY vet Hm, ) r+ 2,04 Y nt Y (A,B) eh otf. (2.5) 


and 


pial 





F,(u,v,r,8) =N,v+ (N,-mx,u,) F+N,b+N,an+N(A,B) +n,+n, (2.6) 


Y(A,B) nordimensional hydrodynamic interaction force 

caused by leading ship. 

N(A,B) nondimensional hydrodynamic interaction moment 

caused by leading ship. 

f,,n, first-order component of random sea. 

f,,n, second-order component of random sea. 

Some assumptions were made for simplification without any 
Significant physical effects. First, the ship’s coordinates 
are assumed to be at the ship’s center of gravity ( Xg, Yo = 
0); second, Y,n and Nn are small terms also neglected from the 
simulation; and third, the second-order force and sway moment 
f, and n, are not considered in this thesis. This is because 
the effect of these slowly varying forces and moments is 
Similar to the interaction forces and moments which are 
analyzed in detail. 

1. Rudder Dynamics Model 

There are two limiters for modeling rudder dynamics, 
the first limiter models the rudder saturation limits at +0.4 
radians. The second limiter models the proportional band of a 
variable-displacement pump by limiting its maximum percent 
stroke, the limits for this nonlinear element have been found 


to be 0.297 radians for this study case. [Ref.1] [Ref.6] 


12 





The lag time and "dead band" in the rudder dynamics 
when a helm angle is commanded are important aspects of the 
manuvering control problem. The rate of change of the rudder 
angle was assumed to be inversely proportional to the error 
signal, subject to a maximum rudder rate of .0698 rad/sec, and 
a maximum error input of 0.297 radians. 

The system gain can be defined as: 

k, = Maximum rudder rate/Maximum error input 

= 0.0689/0.298 = 0.234 /sec 

We nondimensionalize this system gain as: 

k,’=k,*L/u,= 0.234*527.8 ft/(15*1.689) ft/sec =4.89 

Where 

L = ship length= 527.8 ft 
u,;= Ship forward speed= 15 knots= 15*1.689 ft/sec 

This again was taken approximately equal to 5 in this 
study; and 7, = 1/k,’ is defined as the nondimensional rudder 
system time constant. 

The rudder "dead band" was set to + 0.5 deg. An analog 
diagram of the rudder dynamics is shown in Figure 2.4 [Ref.8], 
which includes the response of the rudder system to step 
commands. The equivalent rudder dynamics are described by a 


first order lag, see Figure 2.5, where 


LS 





§=- +644. (2.7) 
na 


T, is the rudder dynamics time constant, and 
6. is the commanded rudder angle. 


2. State Space Representation 


The mathematical model consists of the previous sway 


and yaw equations of motion 


(m-Y,) v+ (mx,-Y,) F=Y,v+ (Y,-mu) r+¥,5+Y(A, B) +f, (2.8) 


(mx,-N,,) v+ (I,-N,) T=N,v+ (N,~mx,u) r+N,5+N(A,B) +n, (2.9) 


and the rudder dynamics 


b=-=2 5+ 28 (2.10) 


The lateral separation rate car be calculated from the 
components of lateral direction of ship forward speed (u) 
lateral speed (v), and heading angle (yp) 


y=usinp+vcosyp (2.11) 
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With the small angle assumption we can linearize this 


equation as: 


y=ytv (2.12) 


where the dimensionless forward speed is one. 

The last kinematic relation is between the ship 
heading y and heading rate per. The first-order wave force f, 
and moment n, in (2.8) and (2.9) vary much rapidly compared to 
the dynamics of the ship and for control/observer design, they 
can be effectively modeled as white noise components. This is, 
however, not true for the interaction forces and moments 
Y(A,B) and N(A,B). These vary at approximately the same rate 
as the ship dynamics since th_, uepend on the relative 
separation between the wo ships. They are best modeled as 
colored noise comp nents: the response of a first-order 


shaping filter driven by white noise 


TyN=-N+Wy (2.13) 


T,Y=-Y+Wy, (2.14) 


where the dimensionless time constants Ty and Ty are equal to 


1(i.e. the time that it takes to travel one ship length). 


15 











Collecting the above equations we get a seven order system in 
state space form. The control input is the commanded rudder 
angle 6. and the augmented state vector contains the actual oe 
rudder angle 6, heading angle ¥, sway velocity v, yaw rate r, 
lateral separation distance y, and the two components of the 
interaction force Y and moment N. In a compact vector 


notation, the complete state space system is written as 


follows. 
X=AX+B,8 .+B.w,+B,w, (22245) 
where 
6 
y 
V 
alae (2.16) 
a4 
N 
Y 
W, 
v.-fee (217) 


16 








(228) 








(2.19) 


[4:3]. 
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(Yy-m) (N,-I,) ~ (Y¥,-mX,) (N,-mX,) 


43,=7 


Y,(N,-I,) -N,(¥,-mX,) 


oe (Y,-m) (N,-I,) - (Y¥,-mX,) (N,-MX¢) 


_ (¥,-mu) (N,-1,) - (¥,-mX,) (N,-mXgu) 


434 (Y,-m) (N,-I,) ~ (Y,-mX,) (N,-mX,) 


Y,-MX_ 


a SS 
36 (¥,-m) (N,-T,) — (¥,-mX,) (N,-mXz) 


N,-I, 


37°" Cy =m) (N,-I,) ~ (¥,-mX_) (N,~mX,) 
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(2.28) 





(2.29) 











__ Y, (N,-mX,) -N, (Y,-m) 











(2.30) 
ad F4 
Y,(N,-MX,) -N, (¥,-m) (2-305 
S6e8 Cr. age 
(N,-mX,) (Y,-mu) ~ (N,-mX,u) (¥-m) (2 32) 
Si. . = oe 
Y,-m (27-33) 
Age = F4 
N,-T, (2332) 
47> F4 
of le (2.35) 
Dae (Y,-m) (N,-J,) -(¥,-mX,) (N,-MmX,) 
iS as es (2.36) 
eT (Y,-m) (N,-I,) -(Y,-mX,) (N-Vv-MX,) 
y,-m (0-37) 
Pua 
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bes (2.38) 


F,=(N,-mX,) (¥,;-mX) - (N,-I,) (Y,-m) (2.39) 


All other terms in the matrices are zero. 
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Figure 2.4 Analog Diagram Representing Rudder Dynamics 


* Note: 
Maximum rudder error signal + 0.29 rad(17 degree). 
K= 0.2/sec. 


Maximum rudder rate 0.0698 rad/sec(3.5 deg/sec). 
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C. EXTERNAL FORCES AND MOMENTS 
External disturbances that act on a ship during underway 
replenishment operations are due to several sources. First, 
the first order wave forces which are rapidly varying forces 
that depend on the seakeeping qualities of the ship and the 
particular seaway during operations. Slowly varying wave 
second order drift forces are important for long term moving 
and positioning operations while the same is true for current 
loads. On the other hand, of particular significance to the 
UNREP problem are the interaction forcesa and moments caused 
by the proximity of the two ships. For this reason we 
concentrate our efforts in the analysis on the control system 
performance of the two categories of forces: First order wave 
and interaction. 
1. Interaction Force and Moment 
These forces arise mainly due to changes in the flow 
field between two ships in proximity to one another and are 
modeled according to the data presented in [Ref.3] [Ref.8]. As 
shown in Figure 2.6 and Figure 2.7 these interaction forces 
are functions of both the longitudinal and lateral separation 
distances of the two ships. Furthermore, they are not of 
consistent signs which makes their effect on path keeping 
control especially troublesome. For simulation purposes the 


data of these two figures are used with bilinear interpolation 
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for the actual values of longitudinal and lateral ship to ship 


separation distances. 


2. First-Order Irregular Wave Effects 


The nondimensional first-order irregular sway force 


and yaw moment are functions of the regular wave encounter 


frequency. The Pierson-Moskowitz energy density spectrum is 


used here aS a computer input to characterize the sea state. 


The nondimensional force and moment can be expressed as [Ref.8] 


¥*(t)=[cos(w,t-e,(w,))/25,(@,) [A,(@,) Fda, — (2-40) 


N*(t)=[cos(w,t-e,(@,)) 2S,(W,) |H,(@,) |*dw, (2.41) 


where 

¥'(t) 

N'(t) 

Q, 

S,(w,) 

E, (w,) 
Em (W,) 

[H, (a) |, | H, (.) | 





is first-order irregular sway force. 
is first-order irregular yaw moment. 
frequency of encounter. 
Pierson Moskowitz energy density spectrum. 
phase angle for sway force. 
phase angle for yaw moment. 
the linear transfer functions at the 
frequency of encounter for the sway force 


excitation and yaw moment excitation. 
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lap co)t2,l aca) are the response amplitude operators (RAO) 


for the sway force and moment at the 


frequency of encounter. 


The data from Table 2.4[Ref.8] were used to determine 


the response amplitude operators and were nondimensionlized 


with respect to the series 60 model. 


For use with the 


Maneuvering model of this thesis a futher nondimensionlization 


step is required. 


For tne nondimensional sway force: 





For the ncudimensional yaw moment: 


where 


M 
H(w,) =¥(@,) —2oe 


1 3 
5 PLsu’ 


L, is the mariner-type ship length. 


L, is the series 60 model length. 


M, is the displacement of series 60 model. 


(2.43) 


The frequency of encounter is function of the absolute 


wave frequency w, the ship speed, the ship-to-wave angle 6, 


and the heading angle y, see Figure 2.8. The equation is 
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e 


w? 
POT. aCoeNe®) (2.42) 


The Pierson-Moskowitz Energy Density Spectrum S,(w) is 
used to characterize the sea state. This is calculated fora 
given significant wave height and for a fully developed, wind- 
driven sea. The characteristics of the wind speed and the 
corresponding significant wave height are shown in Table 2.2. 


A model for Pierson-Moskowitz Spectrum is 


2 )4 
(am) 


2 8 
S,(o) =*Z e (2.43) 
cv) 


where 
a= 8.1 * 10° 
B= -0.74 
g= 32.174 ft/sec’ 
w= frequency of the wave in rad/sec 
V= wind speed in ft/sec 
The corresponding Pierson-Moskowitz Spectrum at the 


frequency of encounter is 


x 


S,(w,) =S, (0) (1-28 ucos (8-w) ) * (2.44) 
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Table 2.2 Characteristics of Series 60 Model 


Length between perpendicular 5.0 ft 1.32 m) 
b 













( 
0.667 ft(_ 0.2m) 
Draft (even keel) 0.267 ft( 0.08 m) 


kg) 


Block coefficient se reece ol) 





Table 2.3 Characteristics of Wind-driven Sea States 


Wind Speed Sea Significant | Average 
(knots) State Wave Height | Period(sec) 








Table 2.4 Nondimensional First-order Sway Force and Yaw Moment in 
Regular Wave 





W, Yaw 
(rad/sec) sway moment 
force phase 
/ (Mg ¢/L} angle 





M= Displaced mass 

g= Acceleration of gravity 

f= Wave amlpitude 

L= Length detween perpendiculars 





Table 2.5 Nondimensional Hydrodynamic Coefficients 


Nondimensional Nondimensional Nondimensional 
Coefficients Factors Values (10°) 





1. The primes ‘' are neglected here but all 
coefficients are dimensionless values. 


2. Xg= 0, L=length between perpendiculars. 


3. u is the ship speed 15 knots, 
nondimensional value is 1. 


4. p is the density of the sea water. 
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FORCE ON SHIP L AS SHIP T PASSES 


Figure 2.6 Ship Proximity Effects: Sway Interaction Forces. 
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igure 2.7 Ship Proximity Effects: Yaw interaction Moments. 
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Figure 2.8 Orientation of Space Axis (%,Y)) and Moving Axis 
(x,y) for Series 60 Model in Oblique Regular Wave 
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III. CONTROL DESIGN 


A. CONTROL METHODOLOGY 

The method for conducting the Underway Replenishemnt at 
Sea maneuver requires the supply ship to maintain a fixed 
speed and heading. The conning officer of the receiving ship 
drives it alongside, and the measurement is the visual 
estimation of the relative distance from the seaman’s eye. 
There are some basic measurement methods and sensor equipments 
which are utilized on the ship and can give very accurate 
distance measurement. The information can be fed into the 
computer which processes the data and automatically maneuvers 
the receiver ship to the desired position. 

A fixed straight line and a zero point were used as the 
course and the center of gravity of the leading ship, and they 
are the references for driving the tracking ship to a position 
of 100 ft from this line, and from -524 ft to 524 ft with 
respect to this point. The controller is at the tracking ship 
(receiving ship). When the tracking ship approaches the 
leading ship (supply ship), the controller of the tracking 
ship will be aware of the excitation due to the interaction 
forces and moments from the leading ship. The dynamics of the 
tracking ship will converge to the sliding surface which is 


chosen in order to drive the ship to the desired position. 
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At the station keeping process, after the tracking ship 


and leading ship are alongside, a constant speed which is a 
normalized speed is considered for both ships. Also constant 
interaction forces and moments are exerted at this stage, 
which will cause a steady state error. Feedforward control is 
used for eliminating this error. 

The feedforward control is mainly used for compensating 
the error of the ship lateral distance at steady state. The 
feedback control law is based on stability requirements, while 
the feedforward part is based on the desired steady state 
accuracy. Since the feedforward control is a function of the 
interaction force and moment, when the latter are zero, the 


former is also zero. Therefore, in general we can write 


Control input(6.)=Feedback Control(6,,) 


+Feedforward Control (8¢,) 


During the control design stage the first-order wave force 
and moment and measurement noise are not considered, and 
perfect and complete state measurements are assumed. Then a 
kalman filter is designed for observing the unmeasurable 
states and the noise filter. 

Figure 3.1 presents a block diagram of the optimal control 


design for the ship replenishment maneuver. 
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B. SLIDING SURFACE DESIGN FOR THE FEEDBACK CONTROL 


1. Reduction of Order 

In sliding mode control, the equivalent system must 
Satisfy not only the n-dimensional state equations but also m- 
dimensional algebraic equations. This indicates that there 
exist only (n-m) indepedent equations on the sliding surface 
and m poles are located at the origin, where n is the number 
of states and m the number of control parameters(1 in our 
case) {Ref.10] [Ref.12]. 


Consider the state equation 


X=AxX+Bu (3.1) 


where states and matrices refer to the mathematical model 
described in Chapter II. 


The first state equation of (3.1) describes the rudder 


dynamics; 


and the rest describe y, v, r, y, Y, and N. 


The system of equations (3.1), (3.2) can be rewritten as: 
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x, _ {An au %|4[P), (3.3) 
X, Az, Az2| |X2} [Bar 
In this form 
x, = [6] 
y 
Vv 
I 
%% = 1yI- (3.4) 
Y 
N 
A, = 0 1 by 6 matrix. 
A,, = Equivalent control 6 by 1 matrix. 
A, = Equivalent state 6 by 6 matrix. 
The linear switching surface parameters are 
K = [K, K,] 
The switching surface can be written as 
p (X) =K,X,+K,x, (3.5) 


where K, =1 and K, is a 1 by 6 matrix which is computed by the 
LQR method for the equivalent system as described later. 


The equivalent control is based on the system 
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%, (t) = [I-Ap, (K,Ap,) Ky) Apo (0) (3.6) 


X,=-Kz Kx, (3.7) 


with x, viewed as the state and x, as the control input. 
By substituting (3.6) into (3.7), the reduced order 


dynamics can be rewritten into the form: 


%,=[4,,-Ai Ky K,) x, (3.8) 


Having found -K,'K, the switching surface and control law is 
completely defined. 

2. ULQR Method for Determination of the Sliding Surface 

The control law design is based on the following: A 5 

degree of rudder for path control is used when the heading 

deviates 5 degrees from zero or the ship reaches a lateral 

offset of 10.43 m (one quarter of beam). From the linear 

quadratic regulator theory the control should minimize the 


quadratic performance function: 


T= = f (x Toxeu TRu) dt (3.95 


where Q is a positive semidefinite matrix weighting on the 


states, and R is a positive definite control weight. 
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The two weights Q and R can be obtained from the criteria 


mentioned above and can be evaluated as follows: 


Weighting on heading angle y: 


( 5 


) -2=131.332 (3.10) 
57.3 





Go2= 


Weighting on side distance y: 


10.43 
290 





Gres | )-25772.463 35) 


Weighting on rudder angle 6: 


2(2 


-2=131.332 3.12 
5703) ( 





Weighting on commanded rudder 6,: 


5 





t= }"*¥21314332 (ac13) 
ie ‘S73 
The equivalent system is: 
X,=A,,X; +A4X, (3.14) 
The Riccati equation is: 
Az3K,+KyAy5~KyAo,0 1A2,K+O=0 (3.15) 
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The equivalent control is then given by: 


X,=-OUAA KX, (3.16) 


3. System Dynamics and Switches 

a. Theory of Switching Control 

Switching control has been successfully applied the 
nonlinear processes such as the "Bang-Bang" control and is 
widely used in optimal control. Systems with switched control 
law represent differential equations with discontinuous right 
hand side. A problem arises in that the conventional 
existence-uniqueness theory for ordinary differential 
equations is not valid. A discontinuous differential equation 
is defined on R" as follows: [Ref.10] [Ref.11] 

Let S be defined as {x:p(x)=0 }, where p is a function 
from R" to R, and S is (n-1) dimensional which represents the 


switching boundary. The switching dynamics are then defined: 


er ce (3.17) 
x=f~ (x) :x=p (x) sO 


where f* and £ are smooth functions from R" to R and are not 
matched on S so that the dynamics are discontinuous at S. The 
mechanism of the switches are operated as p(x) changes sign 
which implies that the state trajectory passes through the 


surface and across it. 
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From the ideal switching law, existence of a sliding 
mode requires infinitely fast switching. But in actual systems 
imperfections exist in the facilities responsible to the 
switching control such as delay, hysteresis, saturation, etc., 
which force the switches with a finite frequency. The ideal 
switching can still be considered as long as the frequency of 
the switching is much higher than the dynamic response of the 
system. This regularization of the switching dynamics is 
schematically depicted Figure 3.2 and Figure 3.3. 

The value p represents the switching variable: when 
pz+l, the dynamics are described by +r, and when ps-i, they 
are described by -r, where r is chosen as 1. A linear relation 
of x and p(x) within this saturation region is described by 
the switching function. 

Existence of the sliding mode requires stability of 
the state trajectory to the sliding surface p(x)=0 and 
asymptotic stablility within the region of attraction which is 


the region between +A and -A. 
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Figure 3.2 Saturation Switch Figure 3.3 Effects of the 
boundary Layer 
Regulation 


b. The Lyapunov Function 
The second method of Lyapunov provides a 
characteristic function for analyzing the existence of the 
sliding mode. The Lyapunov function v(t,x,p) is continuously 
differentiable with respect to all of its variabless and 


satisfies the following conditions: [Reft.12] 


@ v(t,x,p) is positive definite with respect to p. 

@ The total time derivative of v(t,x,p) for the system has 
a negative supremum for all xeQ except for x on the 
Switching surface where the control inputs are undefined. 

The structure of the function vi(t,x,p) is 


determined by the ease which one can compute the switching 


parameters, such as pole placement, Root locus, LOR, or other 





se 





‘ 


state space control theory, for application of variable 
structure control design. For all single input system a 
suitable Lyapunov function is v(t,x) = 0.5p?(x), which is a 
function of p depending on the control and has a negative 


derivative function with respect to time 


dp* _, 2 64 (3.18) 





In the domain of attraction the state trajectory converges to 
the surface and is restricted to the surface for all 
subsequent times. The feedback gains associated with the 
optimal design are computed from Linear Quadratic Regulator 
theory. The weighting matrices for control and state 
associated with the cost function are obtained from (3.10) to 
(35,23 xc 
4. Feedback Control Law 

The method of determination of the switching surface 
p(x)=0 is by Lyapunov function. The stability function 
guarantees that a sliding mode exist at every tet), where t, is 
the time when the state trajectory intercepts the surface. 


The existence of the sliding mode implies tl.at: 


® dp(x)/dt < 0. 


x)=Kx =0, where the K is the surface parameters 


® p( 
(feedback gains), and x is the state vector. 
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From the chain rule 


p (x) = SE X= KK<O (3.19) 


we substitute the constraint equation 


X=AX+Bu (3.20) 


and the definition 


p=-n*sign(p (x) ) (3.21) 


where y is the weighting for the switching structure. 
This form guarantees the existence of a sliding mode and will 


operate as the sign of the state trajectory changes. 


Combining (3.19), (3.20) and (3.21) 
GeI*) x=-n?8ign(p (x) ) (3.22) 


where dp(x)/dx =K and supstituting (3.22), then the control 


law is written as 


— dp (x) -1 dp (x) _ dp (x) =teaidene (3.23 
u a ae Ax- ( ax B) “+n sign(p(x)) , 
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There are two parts in this feedback control. They are 
a feedback structure term for the first part, and the second 
part is for the switching structure. 

In general the steps for the control design are: 

@ Determining the reduced order dynamical equations 
governing the system motion on the switching surface. 

@® Choosing the switching surface parameters K for a linear 

switching surface p(x)=0, so that the system is stake in 


the sliding surface that has been chosen. 


@® Augmenting the control law. 


5. Switching Surface Design 

Generally, the stability of the closed loop system 
depends on the design of the switching surface. The design of 
the switching surface depends on the choice of the surface 
parameters and the dynamics of the system equation. 

The switching surfaces are designed such that the 
state trajectory is restricted to p(x)=0, and that stability 
is guaranteed. 

Again from the reduced order dynamics tne state 


equation can be written as: 


X,=A,, x, +A,,x,+3,,uU 


(3.24) 
X,=A,,%,+A,,X,*B,,u 








The rudder dynamics, first equation of (3.3), will be 


responsible for the control phase of the second equation of 
(3.3), which has an order n-1. The equivalent control gives se 


the equivalent system in the form 


X, =A, 2X tAy X (3.25) 


where 
x, ls the (n-1)th-order equivalent state vector. 
x, is the 1lst-order equivalent control input. 
A, is the equivalent state matrix (n-1) by (n-1) 
matrix. 
A, is the equivalent control matrix (n-m) by m matrix. 


The gain K, results by minimizing the cost function 


J, =[" [x,"0,x,+xR,x,] dt (3.26) 
where 
R,=(@i1] (3.27) 
0,=diag(|qj, 0 0 ges 0 OJ) (3.28) 
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The weighting factors are as indicated by (3.10), (3.11), and 


(3012), (3213). 

Since (A,,A,) is controllable, it is possible to 
effectively use classical feedback control techniques to 
compute a K, and K, such that (A,)+A,,K,'K,) has the desired 
characteristics. Having found [K, K,], the linear switching 


surface is: 


[K, K,]x=0. (3.29) 


C. SLIDING SURFACE DESIGN FOR FEEDFORWARD CONTROL 
1. Introduction 

There are two types of disturbances that were 
considered in this thesis. These are time varying low 
frequency disturbances at passing process, and constant 
disturbances at station keeping process. 

A steady state error in the station keeping process 
indicates that feedback control is not enough. The 
compensation is achieved by feedforward control. 

The design objective is to eliminate the steady state 
error. Prediction of disturbance and cancellation process is 
done by the feedforward control generated at time t. The 
disturbances were assumed to be the same at one time step 


before (t-at) and one time after (t+at) 
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ay=Y(ttat) -Y(t) (3.30) 


AN=N(tt+at) -N(t) 


AY=Y(t)-Y(t-at) (3.31) 


AN=N(t) -N(t-atc) 


By the assumption that at is small and we have slowly 


varying interaction forces and moments as the tracking ship 


passes the leading ship. The differences will approach to zero 


lim aY=0 
at~0 (3.32) 
lim aN=0 


at-0 


2. The Feedforward Control Law 


The steady state Conditions are: 








The separation distance is assumed to be y=0, which is the 


desired position in which the tracking ship will be located. 


We substitute (3.33) into the state equations and we get 


a,5+b,v=c, 


a,5+b,v=c, 


where 


a, =Y, (N,-MX 1 -N, ( ¥=m) 


b,=Y,,(N,-mX,) -N,. (¥,-m) 


C,=-(¥+#,) (N,-mX,) + (¥,-m) (+n) 


a,=Y, (N,-I,) -N, (¥,-mX,) 


b,=Y,.(N,-I,) -N,(¥,-mX,) 


C=7(N AT) (Yet ye Yyomx) (Nen;) 
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.34) 


35) 


.36) 


«37) 


.38) 


39) 


-40) 


-41) 








Simultaneous solutions of (3.34) and (3.35), provide 


the state variables at steady state. The results are: 


_ te, xd) =(c3xb,) 


= 3.42 

S (a,xb,) ~(a,xb,) 
_ (¢,xa,) ~(C,xa,) 

ve" 7b xa = (Deka) (3.43) 

WoH-Vs (3.44) 


The results of these algebraic equations provide the 
amount of the control needed. 
3. Switching Surface Design For Feedforward Control 
As mentioned above rudder angle, heading angle, and 
lateral separation rate are depend on the interaction forces 
and moments at steady state. If constant disturbances act on 
the ship then these three states will not go to zero anda 
steady state error will be presented. An appropriate 
compensation input will be found to eliminate this error. 


The steady state switching surface is: 


P.(X) =K,8.+K,,W.+K,,V,+K,.N+K,.Y (3.45) 








4. Augment of Automatic Control Law 
A subsliding surface was designed for this feedforward 
control input. (An subscript s indicates that the states are 
for feedforward calculation) 


Substitute the feedforward control switching surface 


(3.45) into the control law and augment the reference input as 


8 .=— (KB) “+ KAx,- (KB) “*n*sign(p,) (3.46) 


where 


6, is the feedforward control rudder input. 


p, is the switching surface for feedforward control 


at steady state condition. 


n is the weighting factor for switching control. 


x. (3.47) 


x, is the state vector at steady state. 


The augmented control law for tracking design is completed 


and the schematic block diagram is shown in Figure 3.4. 
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* Note for Fig 3.4: 


B,=(N,-mx,) Y,- (Y,-m) N 


Vv G Vv Vv Vv 


B,=(N,-I,) Y,- (¥,-mx,) N, 


A, =(N,-MX,) Yy- (Y,-M) N, 


A= (N,~I,) Yg- (Y,-mx,) Ns 


f,,n, : 1st-order wave force and moment. 


N,Y : Interaction forces and moments. 
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5. Simulation at Calm Sea 

a. Introduction 

Fortran Program SHIP4QQ.FOR(see Appendix B) was used 
to simulate the UNREP in the horizontal plane using controller 
designed by LQR. Simulation is performed with a perfect and 
full state feedback at this stage. Interaction forces and 
moments are calculated by the bilinear rule. The surge 
velocity is simulated as 15 knots. The time step is .005 
dimensionless seconds which corresponds to 10 Hz sampling 
rate. Two data files TABLE1.DAT and TABLE2.DAT contained the 
data points from Figure 2.6 and Figure 2.7. 

b. Simulation Condition Specification 

Two types of ship maneuvering processes are 
considered in this control test. These are the passing 
maneu*.r and the station keeping process. The passing maneuver 
is a process in which the tracking ship is passing the leading 
ship at constant speed while keeping a desired course. The 
station keeping process is a process in which interaction 
forces and moments are constant after both ships are driven 
alongside, where the same course and speed for the leading and 
tracking ship are maintained at this stage. 
In this chapter, the sea state is considered as a 

calm sea, Perfect measurements from the sensors are assumed. 
The initial conditions for the tracking ship are specified as: 


® One normalized separation distance, y(0)=1. 
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@® Initial heading angle 5 degrees, ¥(0)=.087 rad. 


@® Simulation starts at -4000 ft longitudinal distance, and 
ends at 8000 ft. 


@® The leading ship travels twice as fast as the tracking 
ship during the passing process. 


The weighting factor 7 in the switching structure 
of the feedback control dominates the sensitivity of the 
control law, the stability of the state trajectory in the zone 
of the attraction, and the capability of the design to remain 
stable under external disturbances. Too small a value 
indicates that the system may be slow or even unstable, too 
high a value may cause a chattering problem. The chattering 


problem will not be discussed here[Ref.14] [Ref.12]. 


D. SIMULATION RESULTS AT CALM SEA 
1. Passing Process, y=1, A=0.12 

Results for this case are presented in Figure 3.4. The 
four graphs in this and all similar figures show clockwise 
from the upper left corner the interaction force and moment, 
the actual rudder angle, the tracking ship lateral deviation 
from the desired position( 100 ft), and the ship heading in 
radians. All variables are presented versus the ship-to-ship 
longitudinal separation with 0 corresponding to the two ships 
alongside. It can be seen from the figure that the control 
objectives are met. The effect of the interaction forces is 
evident by the rudder activity at about 0 ft longitudinal 


distance. 
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2. Passing Process, y=1, A=0.06 
The effects of reducing the boundary layer thickness 
are shown in Figure 3.5. The results are virtually identical 
to those of Figure 3.4 which means that for the above range of 
variation, A has no significant effect. It should be 
mentioned, though, that this is true under the current 
assumptions of calm sea and perfect feedback, and it may no 
longer be the case under the general conditions studied in the 
next chapter. 
3. Passing Process, n=8, A=0.12 
The effects of increasing the switching gain 7 are 
shown in Figure 3.6. In order to demonstrate this effect, the 
commanded rudder angle is presented instead of the interaction 
forces. It can be clearly seen that increasing » results in 
more control chattering and increased overshoot of the 
commanded path. Beth of these are operationally dangerous and 
should be avoided in practice. 
4. Passing Process, y=0.5, A=0.12 
The effects of reducing the value of the switching 
gain are shown in Figure 3.7. The control response is now very 
sluggish and the objectives of the underway replenishment 


maneuver are not met. 
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Figure 3.5 UNREP Passing Process, No Feedforward, Calm 
Sea(Prefect Feedback, y=1, A=0.12) 


oF 




















2019 4009 - 4900 -2000 





Figure 3.5 UNREP Passing Process, No Feedforward, Calm 
Sea(Perfect Feedback, n=1, A=0.06) 
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Figure 3.7 UNREP Passing Process, No Feedforward, Calm 
Sea(Perfect Feedback, 7=8, A=0.12) 
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Figure 3.8 UNREP Passing Process, No Feedforward, Calm 
Sea(Perfect Feedback, 94=0.5, A=0.12) 
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5. Station Keeping Process, y=1, A=0.12, No Feedforward 
During the station keeping process where the two ships 
travel alongside for some time, constant interaction force and 
moment will develop, see Figure 3.8. As a result in the 
absence of the feedforward term, a nonzero steady state error 
will be present. Although the feedback control stabilizes the 
system, it cannot ensure the required steady state accuracy. 
6. Passing Process, y=1, A=0.12 
The simulations that follow include the feedforward 
term in the control lawas developed previously. The results 
are shown in Figure 3.9. Comparing these to Figure 3.4. where 
feedforward was not included we can see that for the passing 
process the effects of the feedforward term are very small 
except for the increased rudder action at 0 ft longitudinal 
position, as expected. 
7. Passing Process, 7n=4, A=0.06 
As seem from Figure 3.10 increasing nand decreasing A 
results in more path accuracy and increased control activity. 
The tighter control law in this case helps however during the 
station keeping process in the following two figures. 
8. Station Keeping Process, yn=1, A=0.12 
The effect that the feedforward term has or path 
accuracy is evident by comparing the results of this 
Simulation Figure 3.11, with those at the no feedforward case, 


Figure 3.8. Although the transient response during the 
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Figure 3.9 UNREP Station Keeping Process, No Feedforward, 
Calm Sea(Perfect Feedback, yn=1, A=0.12) 
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Figure 3.10 UNREP Passing Process, Calm Sea(Perfect 
Feedback, Disturbance Compensation, n=1, A=0.12) 


63 








| 





[F 


y 
2003 4000 - 4000 - 2000 0 2000 4008 




















Figure 3.11 UNR&P Passing Process, Calm Sea(Perfect 
Feedback, Disturbance Compensation, 7=4, A=0.06) 
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approach phase is identical, the steady state path deviation 


becomes much smaller and asymptotically reaches zero. 
9. Station Keeping Process, 7=4, A=0.06 
The results of this simulation are shown in Figure 
3.12 where we can clearly see the increased rudder action and 
path overshoot. These result from increasing 7 and decreasing 


4 as in the passing process. 
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Figure 3.12 UNREP Station Keeping Process, Calm Sea(Perfect 
Feedback, Disturbance Compensation, n=1, A=0.12) 
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Figure 3.13 UNREP Station Keeping Process, Calm Sea(Perfect 
Feedback, Disturbance Compensation, 9=4, A=0.06) 
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E. REMARKS 

A sliding mode controller has been designed for the UNREP 
which utilized a LOR method. The controller have been proved 
to be very effective for a time varying disturbance (passing) . 
Disturbance compensation was achieved by means of a 
feedforward term in the automatic controller. 

The results demonstrate excellent tracking and path 
keeping characteristics in the presence of strong interaction 
forces and moments without any compromise in the stability and 
robustness properties at calm seas. 

In this control test results all states are measurable, 
but in reality for a ship operating at sea not all states can 
be measured. So a Kalman filter design is needed as done in 


Chapter IV. 








A. INTRODUCTION 


1. Disturbances 

As was demonstrated in the previous chapter the 
optimal sliding mode control law was able to guarantee 
Stability and steady state accuracy. For that it required 
knowledge of all states including the interaction forces and 
moments. An observer is therefore necessary in order to 
provide an estimate as accurate as possible of all system 
states. The first order wave forces can be effectively modeled 
as white noise components since they vary very rapidly 
compared to the dynamics of the maneuvering ship. A shaping 
filter is introduced in order to model the interaction force 
Y and moment N as exponentially correlated noise driven by 


white noise. In other words 


Ty Ty an 
Yes ye Y 
Ty Ty 


where wy, wy are white noise components and Ty and T, are time 
constants of the filter which are approximately equal to one 


(the time that it takes to travel ine ship length). 
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The power spectral density of Y and N is 


? 0] (4.2) 
vilo 0, 


where Q, and Q, are estimated using 
On= 2 ONT y (4.3) 


= 2027 
Oy=207Ty (4.4) 


The root mean aquare values of N and Y denoted by oy, and 
og, can be estimated using the typical curves for the 


interaction forces shown in Figure 4.1. This calculation gives 


1.548x1078 0 
0.7 (4.5) 


0 8. 97x10" 


2. Measurement Noise 
The measurable states in this study are heading 
angle(y), yaw rate(r), and lateral separation distance(y). 


The output equation is 
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where y, 1S the output vector, the output matrix 





and V, is the measurement noise. The power spectral density R, 


of V, can be evaluated using the numerical values in Table 
4.1. The weighting matrix R, 18 written in the form 
Tv,, 0 0 
R,=|0 rv, 0 (4.7) 
0 0 frv,, 
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Table 4.1 The Power Spectral Density of Measurement 
Noises 





Source 


Measurement 
v gyro 
compass 






IV; = 
4.612*10% 






IVy = 
3.218*10" 







IV33 = 
4502*10" 


o: Random standard deviation frer the true 
value( units are degrees for neading angle 
and yaw rate, ft for y). 

r: The correlation time 
All values of R, are jn nondimensional 
values. 
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B. DESIGN OF KALMAN FILTER 

In order to avoid enlarging the plant noises and 
measuremant noises, an optimal stochastic observer is 
required. The state vector and the covariance error are 


described as follows 


E[x(t,) 1 =X, 


— nee (4.8) 
ELIX(E,) x) Let.) =x, Sp, 
The estimate of the state at ty is assumed to be 
R(t) =X, (4.9) 
The disturbance in the state equations has 
Elw(t)] =0 4.10) 
E[w(t)w7(t)] =O(t)8(t-t) a 
The measurement noise has 


Elv(t)vi7(t))] =R(t)6(t-t) 


Design of the Kalman Filter is based on the cost function 
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Jes [ (x )-Xp) TPo? (x p-X) 1 


+f CwTO-lw+ (V¥_~Cx) TR} (y_-Cx) ] dt 
24 ty 


where the first term minimizes the error in the initial 
estimate, the second term minimizes the effect of w, and the 
third term minimizes the error in the estimate of x. 


The Kalman Filter has the familiar form 


GE -AR+BUtL (Yq~CR) (4.13) 


L is the time-varing Kalman gain matrix, 


L=PCTR-1 (4.14) 


P is the solution of the Riccati differential equation, 


P=AP+PA™+B,0B,’-PC™R CP (4.15) 


B, is the disturbance matrix; and 
Pp (tg) =Po. 
At steady state P is the solution to the algebraic Riccati 


equation 


AP+ Pa™+B,OB,+PC™R*CP=0 (4.16) 
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In such a case the Kalman filter gain matrix is time- 
invariant and the Kalman filter is a full order observer with 


its gains tuned in an optimal way to the system at hand. 
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Figure 4.1 Turning Phase UNREP Passing Process, Open Sea, 
Disturbance Compensation y=1, A=0.12 
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Figure 4.1 Turning Phase- (continued) 
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Figure 4.1 Position Phase(Steady State) 
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Figure 4.1 Position Phase- (continued) 
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Figure 4.2 Turning Phase UNREP Passing Process, Open Sea, 
Disturbance Compensation n=4, A=0.06 
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Figure 4.2 Turning Phase- (continued) 
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Figure 4.2 Position Phase(Steady State) 
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Figure 4.2 Position Phase- (continued) 
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C. SIMULATION AT OPEN SEA 
1. Introduction 
In the following we present results for both the 
passing and the station keeping process. The sliding mode 
linear quadratic Gaussian compensator is used along with the 
feedforward term developed in chapter III. All results are 
based on UNREP in a random sea of significant wave height of 
16 ft. Solid curves in the graphs denote the actual values of 
the variables while dotted curves correspond to the estimated 
values from the Kalman filter. 
2. Passing Process,A=0.12, y=1 
The results of the simulation are shown in Figure 4.1. 
Convergence to the desired track is achieved. The steady state 
deviation of the tracking ship from the commanded path is 
minimal. The steady state use of the rudder angle is also very 


Satisfactory. 


3. Passing Process, 7=4, A=0.06 
The results of the simulation are shown in Figure 4.2. 
Initiai conditions were selected on the desired path so this 
Simulation demonstrates the effects of the passing ship 
disturbance. The two ships are located beam-to-beam at 0 ft 
longitudinal distance. The maximum deviation of the tracking 


ship is + 3 ft. 
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Figure 4.3 Turning Phase Station Keeping Process, Open Sea, 
Disturbance Compensation n=1, A=0.12 
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Figure 4.3 Turning Phase- (continued) 
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Figure 4.3 Position Phase(Steady State) 
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Figure 4.3 Position Phase- (continued) 
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Figure 4.4 Turning Phase UNREP Station Keeping, Open Sea, 
Disturbance Compenastion 7=4, A=0.06 
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Figure 4.4 Tuning Phase- (continued) 
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Figure 4.4 Position Phase(Steady State) 
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Figure 4.4 Position Phase- (continued) 








4. Station Keeping Process, yn=1, A=0.12 


The results of this simulation are shown in Figure 
4.3. Steady state accuracy is maintained. There appears to 
exist an error in the estimates of Y and N despite the zero 
steady state error in the path deviation. This is attributed 
to the existence of the first order wave force and moment. 
Since there is no explicit estimation scheme for these, their 
effect is effectively combined by the observer into an 
equivalent interaction force and moment which of course differ 
from the actual. 

5. Station Keeping process, y=4, A=0.06 

The results of this simulation are shown in Figure 
4.4, Compared to the simulation shown in Figure 4.3 we can see 
the effects of the tighter control law in this case. Both the 
transient and steady state path deviation are reduced while 


the rudder activity is slightly increased. 


D. CONCLUSIONS AND REMARKS 

The sliding mode Linear Quadratic Gaussian compensator was 
able to adequately control the lateral motion of a surface 
ship during underway replenishment in an open sea. The 
feedforward term ensured state path accuracy in the presence 
of interaction forces and moments as well as first order wave 
effects. The Kalman filter was able to minimize the effects of 
the measurement noises and external random disturbances while 


at the same time it provided the controller with an accurate 
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estimate of the states needed for control. Different control 


Options were considered with regards to their transient and 
Steady state characteristics. These comparisons provided 
insight towards appropriate selection of the boundary layer 
thickness and switching gain of the sliding mode control. 
Suggestions for further research include the application 
of integral control instead of disturbance estimation and 
compensation for eliminating the steady state path error. 
Furthermore, the effects of constant disturbances, such as 
Currents, and slowly varying second order wave draft forces 


should be investigated. 
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APPENDIX A. SUMMARY FOR COMPUTER PROGRAMS 


The complete listing of the computer programs written and 
used in this thesis are in APPENDIX B and APPENDIX C. 

The MATLAB program in APPENDIX B is used for calculating 
feedback gains and observer gains using the LQR method, which 
are utilized for sliding surface design and Kalman filter 
design. 

The FORTRAN program in APPENDIX C is a complete listing 
for PASSING PROCESS simulation. The modified program for the 
STATION KEEPING PROCESS is not listed here, the only change is 
the following FORTRAN code logic expression: 


if (b.ge.0) b=0. 
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APPENDIX B. MATLAB LQR DESIGN PROGRAM FOR UNREP SLIDING MODE 
CONTROLLER AND KALMAN GAINS a 


mass=0.0; 
xg=0.0; 

iz=0.0; 

tr=.2; 
nvdot=-19.7e-5; 
xudot=-850e-5; 
xu=-120e-5; 
yv=-1243e-5; 
yvdot=-1500e-5; 
yr=-510e-5; 
yrdot=-27e-5; 
ydr=270e-5; 
nv=-35le-5; 
nr=-227e-S; 
nrdot=-68e-5; 
ndr=-126e-5; 


xn=4.62e-5; 
yne-.52e-5; 
mn=.26e-5; 


apl=yvdot*nrdot -vrdot *«nvdot 

a32=(yrdot*rwv yv*nrdot) /apl; 

a33=(yrdot*nr-yr*nrdot) /apl; 

a35=yrdot/apl; 

a36=-nrdot/apl; 

ap2=nvdot*yrdot -nrdot*yvdot 

a42=(nv*yvdot -nvdot*yv) /ap2 

a43=(nr*yvdot-yr*nvdot) /ap2 

a45=yvdot /ap2 

a46=-nvdot/ap2 

b21=(ndr*yrdot-ydr*nrdot) /apl; 

b31l=(ndr*yvdot-ydr*nvdot) /ap2 

a=[0 0100 0;0 a32 a33 0 a35 a36;0 a42 a43 0 a45 a46; 
110 000;0000-10;0 0000 -1i)j; 

b=[0 b21 b31 0 0 Oj’; 

aab=b31 :a32-b21*a42; 
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aal=(a42*a35-a45*a32) /aab; 
aa2=a35*a42/aab; 

bb1= (a35*b31-a45*b21) /aab; 
bb2=a36*a42/aab; 

tau=.2; 

ad=1; 

all=[-1/tau]; 

al2=[0 0000 0]; 

bb=[{1/tau 00000 0)’; 

aa=[{all al2;b a]; 

gq=diag([131.3 0 0 771.69 0 0]); 
qq=diag([131.3 131.3 0 0 771.69)); 
r=[(131.3]; 
c='100000;001000;00010 0); 
{k,s]=lqr(a,b,q,r); 

ss=(1 k] 

gg=inv(ss*bb) *ss*aa 
kn=inv(ss*bb) ; 

ca=aa-bb* (gg-kn1) ; 

cc=[0 100000;0001000 
gl1=[0 000010;000000 
qi=[1.548e-8 0;0 8.97e-8]; 
rr=diag([{4.612e-8 3.213e-7 4.502e-6]); 
l=lqe(aa,gl,cc,qi,rr) 

f=aa-1l*cc*taa 

bf=bb-1*cc*bb 


7000010 0]; 
a) * 3 
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APPENDIX C. UNREP SIMULATION PROGRAM: SLIDING MODE 
CONTROLLER DESIGNED BY LOR METHOD 


program shipdynamics 
CORR KKK KK KK KKK KK KKK KKK KKK KKK KKK KKK KKK KKK KKK KEKKEKEKKKEKEKKEKKKKKEK 


Simulation program for UNREP 
Written by LT FU, HSU-SHENG SIMON 


The program controls one rudder input with a sliding mode 
design for path keeping and tracking. 


Speed is constant 1 normalized value. 


The initial conditions for ship are follows: 1 
nondimensional leteral seperation distance, -4000 ft 
longitudinal distance, 0.087 rad heading angle. 


The disturbances include interaction forces and moments, 
and first-order wave sway force and yaw moment. 


The observed states: yaw rate, lateral separation rate, 
interaction forces and moments. 


The measurements: heading angle, lateral separation 
distance, and rudder angle. 


The outputs of the program are written to files, these 
files are then accessed by the MATLAB to generat output 
graphics. 


This program utilized subroutines for different functions 
that are required: 
subroutine locate: find the interaction forces and 
moments at the real time and use 
subroutine trlint as interpolation. 
subroutine trlint: incooprate with subroutine locate. 


gQaqgqaqaaqaagnaagaannnqaaagaanaaaaaqaaaqa¢aadqraanaaqaaaaaaaa 
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subroutine randm : generate random signal. 
subroutine trap : Trapzoid integration. 
Subroutine first : first order wave calculation. 


Two data files for interaction forces and moments: 
Tablel.dat: interaction forces. 


Table2.dat: interaction moments. 
CEH KKK KKK KEKE KKK KKK KEKE KKK KKK KEKE KKKKKKEKEKEKKKKKKKKKKEK 


qQgqaqaqaaa 


real mass,nrdot,nvdot,nr,nv,ndr,nn 

real k1,k2,k3,k4,k5,k6,k7,kn,iz,mo,a,b,1 

real w,coef,f1,n1,wn,wy,we2(5),ws 

real o1(7,3),ae(7,7),ba(7,1),bc(7,2) ,bd (7,2) ,sf2(5), 
1 ym2 (5) 


real time 


c 
dimension aa(23) ,bb(11),table1(23,11) ,table2 (23,11) 
data aa/-550,-500, -450, -400, -350, -300, -250,-200,-150, 
1 -100,-50,0,50,100,150,200,250,300, 350,400, 
1 450,500,550/ 
data bb/50,60,70,80,90,100,110,120,130,140,150/ 
data 01/0,2.6383,-1.4846,6.6537, .165, .0167, -.0133, 
1 0, .9536,-3.8899,21.5836,-.911, .1926, -.1238, 
1 0, .0017,3.3203,-.0651,2.566, .6107, .0507/ 
c 
m =23 
n =11 
1 =527.8 
iseed =17 
c 
mass = 0.0 
xg = 0.0 
iz = 0.0 
mvdot=- 19.7*1.e-5 
tr = .2 
xudot=- 850.0*1.e-5 
xu =- 120.0*1.e-5 
yv- =-1243.0*1.e-5 
yvdot=-1500.0*1.e-5 
yr =- 510.0*1.e-5 
yrdot=- 27.0*1.e-5 
ydr = 270.0*1.e-5 


99 











Q 


70 
60 


nv =- 351.0*1.e-5 
nr =- 227.0*1.e-5 
mrdot=- 68.0*1.e-5 
ndr =- 126.0*1.e-5 
xn 7 4.62*1.e-5 
yn =- S52*1.e-5 
nn 7 26*1.e-5 


define system matrix 


do 60 i <=1,7 

do 70 j =1,7 

ae(i,j)=0.0 

ba(i,1)=0.0 

be (i,1)=0.0 

bd(i,1)=0.0 

be (i,2)=0.0 

bd (i,2)=0.0 

continue 
continue 
ael = nrdot*yvdot-yrdot*nvdot 
ae2 = nvdot*yrdot-nrdot*yvdot 
ae(1,1)=-1/tr 
ae(2,4)= 1. 
ae(3,1)= (ndr*yrdot-ydr*nrdot) /ael 
ae(3,3)= (yrdot*nv-yv*nrdot) /ael 
ae(3,4)= (yrdot*nr-yr*nrdot) /ael 


ae(3,6)= yrdot/ael 

ae (3,7) =-nrdot/ael 

ae(4,1)= (nar*yvdot-ydr*nvdot) /ae2 
ae(4,3)= (yvdot*nv-yv*nvdot) /ae2 
ae(4,4)= (nr*yvdot-nvdot*yr) /ae2 
ae(4,6)= yvdot/ae2 

ae (4,7) =-nvdot/ae2 


ae(5,2)= 1. 
ae(5,3)= 1. 
ae(6,6)=-1. 
ae(7,7)=-1. 
ba(i,1)=-1/tr 
be(6,1)= 1. 
be(7,2)= 1. 


100 


d(3,1)= Basra Hage 
d(3,2)=-nrdot/ael 
ae 1)= Baoe shes 
bd (4,2) =-nvdot/ae2 


open(10, file=’sim.dat’,status='old’ ) 
open(11,file=’siml.res’,status='’new’ ) 
open(12,file='’sim2.res’',status=’ new’ ) 
open (13,file='’sim3.res’,status=’ new’ ) 
open (14,file=’tablel.dat’,status='old’) 
open(15,file=’table2.dat’,status=’old’) 
open(16,file=’sim4.res’,status=’new’ ) 
open(18,file='’sim5.res’,status=’ new’ ) 
open(19,file='’sim6.res’,status='new’ ) 
open (20,file=’simol.res’,status=’new’ ) 
open(21,file='’simo2.res’,status=’'new’ ) 


c100 
do 40 i=1,m 
read (14,30) (tablel(i,j),j=1,n) 
read (15,30) (table2(i,j),j=1,n) 
30 format (11£7.2) 
40 continue 


read(10,*)stime,delta 

read (10,*)s1,s2,s3,s4,s5,s6,8s7 
read(10,*)k1,k2,k3,k4,k5,k6,k7,kn 
read(10,*)sigsat,eta,etal,psi,y,xref,yref,ws 


wS=ws* .51444/.3048 


Q 


define state v,r,dr,drc,psi,y 


isim=stime/delta 
v =0.0 


QO, 

K 

W 

[o) 
ooo0 © 
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simulation start 


do 1 i=1,isim 
if (b.ge. (3.5*xref)) stop 


First order wave 
call first (psi,ws,time,f1,n1) 
White noise 


call randm(iseed,wn) 
call randm(iseed, wy) 
wn= (wn*1.e-5) 
wy= (wy*1.e-5) 


Interaction Force and Moment 


a =yref+y*l 

b =x*l1-xref 

call locate (bb,n,a,nx) 

call loctae(aa,m,b,my) 

call trlint(m,n,my,nx,tablel,bb,aa,a,b,yf) 
call trlint(m,m,my,nx,table2,bb,aa,a,b,mo) 
yf=yf*1.e-5 
mo=mo*l.e-5 

rE ble t= 550). 6r-bigt<550) then 
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yf=0.0 
mo=0.0 
endif 


al=mass-yvdot 

bi=mass*xg-yrdot 

Cl=yv*v+ (yr-mass) *r+ydr*dr+yf+f1 
a2=mass*xg-nvdot 

b2=iz-nrdot 

c2=nv*v+ (nr-mass*xg) *r+ndr*dr+mo+nl 
a3=mass-xudot 


Ardot =-dr/tr+drc/tr 


psidot= r 

vdot = (cl1*b2-c2*b1) /(al*b2-a2*b1) 
rdot = (cl*a2-c2*al) /(a2*bl-al*b2) 
ydot = sin(psi) +v*cos (psi) 

xdot = cos(psi) -v*sin(psi) 


yfdot =-yf+wy 
modot =-mo+wn 


Euler Integration for Actual States 


dr =dr +drdot *delta 
psi=psi+psidot*delta 
v =v +vdot *delta 
r =r +rdot *delta 
y =y +ydot *delta 
x =x +¢+xdot *delta 
yf =yf +yfdot *delta 
mo =mo +modot *delta 


Measurement Noises 


call randm(iseed,rnpsi) 
call randm(iseed,rnr) 
call randm(iseed, rny) 
rnpsi=rnpsi*0.0035 

rnr -=rnr *0.0002 

rny =rny *0.0063 
psime=psi +rnpsi 
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a 


a 


1 


rme =r +rnr 
yme =y +rny 


Feedforward Control Rudder Angle Design - 


aal = ydr*nvdot -ndr*yvdot 

bbl = nvdot*yv-yvdot*nv = 
ecd =- (oy£+f£1) *nvdot+ (omo+n1) *yvdot 

aa2 = ydr*nrdot-ndr*yrdot 

bb2 = yv*nrdot-nv*yvdot 

cc2 = (omo+nl) *yrdot - (oy£+f1) *nrdot 

arc0O = (ccl*bb2-cc2*bb1) / (aal*bb2-aa2*bb1) 

v0 = (ccl*aa2-cc2*aal) / (bb1*aa2-bb2*aal) 

psiO =-v0 


sigmal=si*drc0+s2*psi0+s3*v0+s6*omo+s7*oy£f 
if (abs(sigmal).1lt.sigsat) satsgni=sigmal/sigsat 


if (sigmal.le.-sigsat) satsgnl=-1.0 

if (sigmal.gt. sigsat) satsgnl= 1.0 

ad =-etal**2*kn*satsgnil 

dro = Grc0+(k1*drc0+k2*psi0+k3*v0+k6*omo+k7*oyf) -a0 


rudder input calculation 


sigma=s1*odr+ (s2*opsi+s3*ov+s4*or+s5 *oy+sé6*omo+s7*oyf ) 
if (abs(sigma).lt.sigsat) satsgn=sigma/sigsat 

if (sigma.le.-sigsat) satsgn=-1.0 

if (sigma.ge. sigsat) satsgn= 1.0 


drchat=- (k1*odr+k2*opsi+k3*ov+k4*or+k5 *oy+k6*omo+ 


k7*oyf) 
drcbar= eta**2*kn*satsgn 
dre = drchat-drcbar+dr0 
drf = drchat-drcbar 


if (drce.ge. 0.4) drce= 0.4 
if (dre.le.-0.4) drc=-0.4 


Observed States 


obl=ae (1,1) *odr+ae(1,2) *opsitae (1,3) *tov+ae (1,4) *or+ 
ae(1,5) *oy+ae (1,6) *omo+ae (1,7) *toyft ‘ 
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cc 


ob2=ae (2,1) todr+ae (2,2) *opsi+ae (2,3) *ov+ae (2,4) *or+ 
ae (2,5) toy+ae (2,6) *omo+ae (2,7) *oyf 
ob3=ae (3,1) todr+ae (3,2) *opsitae (3,3) tov+ae (3,4) *or+ 
ae (3,5) *oy+ae (3,6) *omo+ae (3,7) *oyf 
ob4=ae (4,1) *odr+ae (4,2) *opsi+ae (4,3) tov+ae (4,4) *or+ 
ae (4,5) *oy+ae (4,6) *omo+ae (4,7) *oyt 
ob5=ae (5,1) todr+ae (5,2) *opsitae (5,3) *ov+ae (5,4) *or+ 
ae (5,5) *oy+ae (5,6) *omo+ae (5,7) *oyf 
ob6=ae (6,1) *odr+ae (6,2) *opsi+ae (6,3) *ov+ae (6,4) *or+ 
ae (6,5) *oyt+ae (6,6) *omo+ae (6,7) *oyf 
ob7=ae (7,1) *odr+ae(7,2) *opsi+ae (7,3) *ov+ae (7,4) *or+ 
ae(7,5) *oy+ae(7,6) *omo+ae (7,7) toyf 
ob23 =bd (3,1) *(n1)+bd(3,2) *f1 
ob24 =bd(4,1)*(n1)+bd(4,2)*f1 
odrdot =obl+ol (1,1) * (psime-opsi) +01 (1,2) * (rme-or) + 
o1 (1,3) * (yme-oy) +ba(1,1) *dre 
opsidot=ob2+o0l (2,1) * (psime-opsi) +01 (2,2) * (rme-or) + 
01 (2,3) * (yme- oy) 
ovdot -=0b3+cl1(3,1)* (psime-opsi) +01 (3,2) * (rme-or) + 
01 (3,3) * (yme- oy) +0b23 
ordot =o0b4+01(4,1)*(psime-opsi) +01 (4,2) * (rme-or) + 
01 (4,3) * (yme- oy) +0b24 
oydot ~ob5401 (5 1) * (psime-opsi) +01 (5,2) * (rme-or) + 
L(53)% gage oy) 
omodot Dee ae 1) * (psime-opsi) +01 (6,2) * (rme-or) + 
o1(6, 38 gma oy) 
oyfdot =ob7+01 (7,1) * (psime-opsi) +01 (7,2) * (rme-or) + 
COLT, a) * tyme oy) 


Eular Integration for Estimate States 


odr =odr +odrdot *delta 
opsi=opsi+opsidot*delta 
ov =Ov +ovdot *delta 
or =or +ordot *delta 
oy =oy +oydot *delta 
oyf =oyf +oyfdot *delta 
omo =omo +omodot *delta 


time=i*delta 
if (mod(i,20).eq.0) then 
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write(11,20)time,psi,v,y,r,dr 


20 format (f10.4,1x,f10.4,1x,f£10.4,1x,f10.4,1x,f10.4, 
1 ix, f10.4) 
write(12,*)x,y,a,b ee 
write (13,*)drchat,drcbar,drc,drf 
write (16,*)yf,mo,dr0 
write (19,*)wy,wn x= 
write (20,50) odr,opsi,ov,or,oy 
50 format (£10.4,1x,£10.4,1x,f£10.4,1x,f£10.4,1x,£10.4) 
write (21,*)omo,oyf 
write(18,*)f1,n1l1 
endif 
1 continue 
stop 
end 
c303 
eccecececeececoceccceececececceecceccceccecctccecececccceccac 
c 
subroutine locate (xx,n,x,j) 
c 
eC This is a subroutine for finding the interaction forces 
c and moments. For location the interaction forces and 
C moments in the input file TABLE1.DAT and TABLE2.DAT 
C correspond the lateral distance and longitudinal 
c distance. 
c 
dimension xx(n) 
c 
do 20 i=l,n 
1f (x.eq.xx(i)) jai 
20 continue 
do 10 i=1,n-1 
1f (x.gt.xx(i).and.x.1lt.xx(i+1)) then 
72h 
endif 
10 continue 
return 
end af 
c 
ccceccecccecececeecceccecccecccceccecceccecececccecececeecece ¥ 
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subroutine trlint(m,n,i,j,ya,xla,x2a,x1,x2,y) 


This subroutine is used for calculating the interaction 
forces and moments, the data get from SUB LOCATE and the 
values calculated by interpolation. 


qaaaaana 


dimension ya(m,n),xla(n) ,x2a(m) 
real t,u,tl,ul 


yi=ya (i 
y2sya (i 
y3= soe 

eds j+1 


t=(x1-xla(j))/(xla(j+1)-x1la(j)) 
u=(x2-x2a(i))/(x2a(i+1)-x2a(i)) 


tl=yl+u* (y3-yl) 
ul=y3+u* (y4-y2) 
y =tl+t* (ul-t1) 


return 

end 
eacccceccoecoccecetacecececetcocceccecececcecceccecececedecece 

subroutine randm(iseed, rndm) 


e 
c generate a gaussian random number 
c iseed - input a integer variable seed 
a output seed is changed during execution for next 
c call 
c rndm - gaussian random variable with zero mean and unit 
c variance 
c 
real rndm,rndl1,rnd2 
: integer iseed,jr,kr,mr 
parameter (jr=5243,kr=55397,mr=262139) 
i iseed=mod (iseed*jr+kr,mr) 
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rndi =(iseed+.5)/mr 

iseed=mod (iseed*jr+kr,mr) 

rnd2 =(iseed+.5)/mr 

rndm =sqrt (-2.*log(rnd1) ) *cos(6.2831853*rnd2) /2.5 
return 
end 


cececececececccececccececececceccececcececccceccceccecccecccecececcccceccceccc 


aQ 


subroutine trap(n*,vx,vy,out) 
numerical integration routine using the trapezoidal rule 


dimension vx(1),vy(1) 
real out 


nl=nt-1 

out=0.0 

do 1 i=1,n1 
outil=.5* (vx (i) +vx(i+1)) * (vy (i+1) -vy(i)) 
out =oOut+outl 

continue 

return 

end 


ceccecceccecccececcccececceccceccccccccccccececccccececcccecce 


9aagaaqaqaaaqnteaaqaeqnqnaeaqaaregaana 


a 


subroutine first (psi,ws,time,f1,n1) 


This subroutine is used for calculating the sway force 
and yaw moment component from first-order wave 
excitation. 


The inputs of the subroutine are: 
(1) The encounter frequency (w,) : 


w.=w- (w~2) /g*ut*cos (theta-psi) 
Ww is the actual wave frequency. 
theta(@) is the wave-to-ship angle. 
psi(y) is the ship heading. 
g is the acceleration of gravity. 
u is the ship forward speed(15 knots). 


(2) The wind speed(ws). 
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(3) The RAO of sway force and RAO of yaw moment are the 
function of encounter frequency 
(4) The Pierson-Moskowitz spectrum is the function of we 
and wind speed 
The associated force and moment are calculated from 
Trapizoid Intergration and be the outputs of the 
subroutine. 
The phase angles correspond to the sway force and 
moment are generated by random process (irregular wave) 
(epsf,epsm) from 0 to 27 


Cc 
Cc 
c 
Cc 
c 
Cc 
c 
Cc 
Cc 
c 
Cc 


c397 

real sf(19),ym(19) ,wee(19) ,w£ (19) ,wm(19) ,vecx(19), 

1 vecm(19) 

real we,alpha,beta,ws,g,epsf,epsm,sx(19) ,psi,w1(19), 
al w2 (19) 

real raom(19),raof(19),time,model,1,1m,rho,ul1,omega, 
1 w(19) ,s (19) 

reai f1,nl 

data sf/1.3191,1.9473,2.5082,2.9529,3.2653,3.3462, 
1 3.0542,2.3640,1.5397, .6624, .21465, .73203,1.033, 
1 1.0222, .68012, .27829, .37778, .55739, .46259/ 
data ym/.07246, .03741, .08306, .22726, .42925, .6717, 

1 .89061, .99661,1.0366, .9459, .73952, .47369, .19894, 

.17588, .31972, .34974, .25216, .07516, .15072/ 

data wee/.293,.361, .433, .509, .588, .679, .756, .845, 

1 -938,1.304,1.133,1.236,1.342,1.452,1.565,1.682, 
1 1.801,1.925,2.052/ 

data wf/274,277,277,280,287,290,293,294,292,281,16€3, 

1 122,110,101, 87,40, 313,290,280/ 

data wm/169,138,32,19,22,23,23,26,27,24,19,10, 339, 

1 242,211,197,185,146,24/ 

c 

model= 33.27 

im ee oe 

1 = 527.8 

rho = 1.9905 

g = 32.2 

ul = 26.0 

alpha= 8.1*.001 

beta =- .74 
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ised = 17 
omega= 2.0 
time = time*l/ul 


c Pierson-Moskowitz Spectrum 
do1i1 ~=#=1,19 -. 
raof (i)=sf (i) *model*g/1m 


i) 
raom(i)=ym(i) *model*g 
wi (i) =(1+sqrt (1-4*ul*cos (2.618-psi)j/g))/ 


1 (2*ul*cos (2.618-psi) /g) 
w2 (i) =(1-sqrct (1-4*ul*cos (2.618-psi) /g) )/ 
1 (2*ul*cos (2.618-psi) /g) 


1f (wl(i)*w2(i).1t.0) then 
lf (wl(i).gt.0) w(i)=wi(i) 
Lf (w2(i).gt.0) w(i)=w2(i) 
endif 
4£ (w1(i)*w2(i).gt.0) then 
if (abs(wl(i)).gt.abs(w2(i))) w(i)=abs(w1(i)) 
if (abs(w2(i)).gt.abs(w1(i))) w(i)=abs (w2 (i) ) 


endif 
s (i) = (alpha* (g**2) /(w(i) **5))* 
1 exp (beta* (g/ (ws*w(i))) **4) 
sx(i) =s(i)/(1- (2*w(i)/g) *ul*cos(2.618-psi) ) 
vecx (1) =sqrt (2*sx(i)*raof(i))}* 
1 cos (wee (i) *time-wf (i) *3.14/180) 
vecm(i)=sqrt (2*sx (i) *raom(i))* 
1 cos (wee (i) *time-wm(i) *3.14/180) 
continue 


Force and Moment calculation 


qaaake 


call trap(19,vecx, wee, f1) 
call trap(19,vecm,wee,n1) 
fl=f1/(.5*rho*tl**2*ul**2) 
nl=nl/(.5*rho*¥l**3*u1l**2) 
return 

end 





as 
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